# -*- coding: UTF-8 -*-
from abaqus import *
from abaqusConstants import *
import regionToolset

session.viewports['Viewport: 1'].setValues(displayedObject=None)

#----------------------------------------------------------------------------------------------
# Create the model
mdb.models.changeKey(fromName='Model-1', toName='Plate Bending Model')
plateModel = mdb.models['Plate Bending Model']

#-------------------------------------------------------------------------------
# Create the part

import sketch
import part 

# a) Sketch the plate using the rectangle tool
plateProfileSketch = plateModel.ConstrainedSketch(name='Plate Sketch', sheetSize = 20)

plateProfileSketch.rectangle(point1 = (0.0,0.0), point2 = (5.0,3.0))

# b) create a shell named "Plate" using the sketch

platePart = plateModel.Part(name= 'Plate', dimensionality = THREE_D, type=DEFORMABLE_BODY)

platePart.BaseShell(sketch = plateProfileSketch)

#--------------------------------------------------------------------------------------------
# Create material 

import material

# Create naterial AISI 1005 Steel by assigning mass density, youngs modulus and poissons ratio

plateMaterial = plateModel.Material(name='AISI 1005 Steel')
plateMaterial.Density(table=((7872, ),  ))
plateMaterial.Elastic(table=((200E9, 0.29), ))

#--------------------------------------------------------------------------------------------------
# Create homogeneous shell section of thickness 0.1m and assign the plate to it

import section

# Create a section to assign to the plate
plateSection = plateModel.HomogeneousShellSection(name='Plate Section',
        material='AISI 1005 Steel', thicknessType=UNIFORM, thickness=0.1)

#assign the plate to this section 
plate_face_point = (2.5, 1.5, 0.0)
plate_face = platePart.faces.findAt((plate_face_point,)) 
plate_region = (plate_face,)

platePart.SectionAssignment(region=plate_region, sectionName='Plate Section',
            offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='')

#--------------------------------------------------------------------------------
# Create the assembly 

import assembly

# Create the part instance
plateAssembly = plateModel.rootAssembly
plateInstance = plateAssembly.Instance(name='Plate Instance', part=platePart,
                        dependent=ON)

#------------------------------------------------------------------------------------
# Create the step 
import step

# Create a static general step
plateModel.StaticStep(name='Load Step', previous='Initial',
        description= 'Apply concentrated fores in this step', nlgeom=ON)

#----------------------------------------------------------------------------------------
# Create the field output request

# change the name of field output request 'F-Output-1' to 'Output Stresses and Displacements'
plateModel.fieldOutputRequests.changeKey(fromName='F-Output-1',
            toName='Output Stresses and Displacements')

# since F-Output-1 is applied at the 'Apply Load' step by default, 'Selected Field Outputs' will be too
# We only need to set the required variables
plateModel.fieldOutputRequests['Output Stresses and Displacements'].setValues(variables=('S','UT'))

#----------------------------------------------------------------------------------------------------
# Create the history output request

# we don't want any history outputs so lets delete the existing one 'H-output-1' del plateModel.historyoutputRequests['H-Output-1']

#--------------------------------------------------------------------------------------------
# Apply boundary conditions-fix one edge

fixed_edge = plateInstance.edges.findAt(((0.0, 1.5, 0.0), ))
fixed_edge_region=regionToolset.Region(edges=fixed_edge)

plateModel.DisplacementBC(name='FixEdge', createStepName='Initial',
        region=fixed_edge_region, u1=SET, u2=SET, u3=SET, ur1=SET, ur2=SET, ur3=SET,
        amplitude=UNSET, distributionType=UNIFORM, fieldName='', localCsys=None)

# Instead of using the displacenents/rotations boundary condition and setting all six DOF to zero
# We could have just used the Encastre condition with the following statement 
# plateModel.EncastreBC(name='Encaster edge', createStepNamem = 'Initial',
        # region=fixed_edge_region)

#------------------------------------------------------------------------------------------------
# Create vertices on which to apply concentrated forces by partitioning part

# create the datum points
platePart.DatumPointByCoordinate(coords=(0.0, 1.0, 0.0)) 
platePart.DatumPointByCoordinate(coords=(0.0, 2.0, 0.0)) 
platePart.DatumPointByCoordinate(coords=(5.0, 1.0, 0.0)) 
platePart.DatumPointByCoordinate(coords=(5.0, 2.0, 0.0)) 

# Assign the datum points to variables
# Abaqus stores the 4 datum points in platePart.datums
# Since their keys may or may not start at zero, put the keys in an array sorted in ascending order
platePart_datums_keys = platePart.datums.keys() 
platePart_datums_keys.sort()
plate_datum_point_1 = platePart.datums[platePart_datums_keys[0]]
plate_datum_point_2 = platePart.datums[platePart_datums_keys[1]] 
plate_datum_point_3 = platePart.datums[platePart_datums_keys[2]]
plate_datum_point_4 = platePart.datums[platePart_datums_keys[3]] 

# Select the entire face and partition it using two points
partition_face_pt =(2.5, 1.5, 0.0)
partition_face = platePart.faces.findAt(partition_face_pt,)
platePart.PartitionFaceByShortestPath(point1= plate_datum_point_1, point2=plate_datum_point_3,
                                faces=partition_face)
                                
# Now two faces exist, select the one that needs to be partitioned
partition_face_pt =(2.5, 2.0, 0.0)
partition_face = platePart.faces.findAt((partition_face_pt,))
platePart.PartitionFaceByShortestPath(point1=plate_datum_point_2, point2=plate_datum_point_4,
                        faces=partition_face)
                        
# Since the partitions have been created, vertices can be extracted
vertices_for_concentrated_force = plateInstance.vertices. findAt(((5.0, 1.0, 0.0),),
                ((5.0,2.0,0.0),),)

# ----------------------------------------------------------------------------------------
# Apply concentrated forces

plateModel.ConcentratedForce(name='Concentrated Forces',
        createStepName='Load Step',
        region=(vertices_for_concentrated_force,),
        cf3=7000.0, distributionType=UNIFORM)
        
#-------------------------------------------------------------------------------------------
# Create the mesh

import mesh

#set element type
plate_mesh_region  = plate_region

elemType1 = mesh.ElemType(elemCode=S8R, elemLibrary=STANDARD)

platePart.setElementType(regions=plate_mesh_region, elemTypes=(elemType1,))

# seed edges by number
mesh_edges_vertical = platePart.edges.findAt(((0.0, 0.5, 0.0), ),
                                            ((0.0,1.5,0.0),),
                                            ((0.0,2.5,0.0),),
                                            ((5.0,0.5,0.0),),
                                            ((5.0,1.5, 0.0),),
                                            ((5.0,2.5,0.0),))
mesh_edges_horizontal = platePart.edges.findAt(((2.5, 0.0, 0.0),),
                                                ((2.5,1.0,0.0),),
                                                ((2.5,2.0,0.0),),
                                                ((2.5,3.0,0.0),))
platePart.seedEdgeByNumber(edges=mesh_edges_vertical, number = 3)
platePart.seedEdgeByNumber(edges=mesh_edges_horizontal, number=10)

platePart.generateMesh()

#-----------------------------------------------------------------------------------
# Create and run the job

import job
# create the job

mdb.Job(name='PlateJob',model='Plate Bending Model',type=ANALYSIS,
            description='Job simulates the bending of a plate')
# run the job
mdb.jobs['PlateJob'].submit(consistencyChecking=OFF)

# do not return control till job is finished running
mdb.jobs['PlateJob'].waitForCompletion()

##################################################################
#---------------------------------------------------------------------------
# Post Processing
#---------------------------------------------------------------------------
##################################################################

#-----------------------------------------------------------------------
# Display deformed state with contours

import visualization

plate_viewport = session.Viewport(name='Plate Results Viewport')
plate_Odb_Path = 'PlateJob.odb'
an_odb_object = session.openOdb(name=plate_Odb_Path)
plate_viewport.setValues(displayedObject=an_odb_object)
plate_viewport.odbDisplay.display.setValues(plotState=(CONTOURS_ON_DEF, ))

# Report stresses in descending order

import odbAccess
# The main session viewport must be set to the odb object using the following line.
# If not you might receive an error message that states
# "There are no active entities. No report has been generated."
session.viewports['Viewport: 1'].setValues(displayedObject=an_odb_object)

# Set the option to display the reported quantity (in our case the stresses) in
# descending order
session.fieldReportOptions.setValues(sort=DESCENDING)

# Name the report and give it a path. If you do not assign a path (as is done here)
# it will be stored in the default abaqus temporary directory
report_name_and_path='PlateStresses.rpt'

#........................................................................................
# You may enter an entire path if you wish to have the report stored in a
# particular location.
# One way to do it is using the following syntax.
# report_name-'PlateReport
# report_path='C: /MyNewFolder/'
# report_name_and_path = report_path + report_name +',rpt'
# Alternatively you could have used 1 statement instead of these 3 :
# report_name_and_path='C:/MyNewFolder/PlateReport.rpt
# Note however that the folder 'MyNewFolder' must exist otherwise you will likely get the following error
#＂IOError: C/MyNewFolder: Directory not found＂
# You must either create the folder in Windows before running the script
# Or if you wish to create it using Python commands you must use the os.makedir()
# or os.makedirs()function
# os.makedirs() is preferable because you can create multiple nested directories in one statent if you wish
# Note that this function returns an exception if the directory already exists
# hence it is a good idea to use a try block

#try:
#   os.makedirs(report_path)
#except:
#   print ＂Directory exists hence no need to recreate it. \
#                                                Move on to next statement＂
# Here it is rewritten without all the comments
"""
report_name='PlateStresses
report_path='C: /MyNewFolder/'
report_name_and_path  report_path ＋ report_name ＋ '.rpt'
try:
    os.makedirs(report_path)
except:
    print ＂Directory exists hence no need to recreate it. \
                                                    Move on to next statement"
"""
#.................................................................................

# Write the field report outputting the Mises stresses
session.writeFieldReport(fileName=report_name_and_path, append=OFF,
            sortItem='S.Mises', odb=an_odb_object, step=0, frame=1,
            outputPosition=INTEGRATION_POINT,
            variable=(('S', INTEGRATION_POINT, ((INVARIANT,'Mises'), )),))